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Abstract 



^0 ' In this report, we address differential systems with Lipschitz non linearities; this 

f] , study is motivated by the subject of vibrations of structures with unilateral springs or 

non linear stress-strain law close to the linear case. We consider existence and solution 
with fixed point methods; this method is constructive and provides a numerical algo- 
rithm which is under study. We describe the method for a static case example and we 
address periodic solutions of differential systems arising in the vibration of structures. 



1 Introduction 

m ; 

This work stems from non destructive testing experiments reported in [VL04J and from the 
experiments and the computations of a beam model of the dynamic of satellite solar arrayqj 
[HFR09j,|Hazl0j. In these two situations the structure is submitted to an harmonic force; 
the behavior of the structure is locally non linear in particular due to unilateral elastic 
contacts. 

During these studies, it became clear that an in-depth understanding of the forced 
dynamic behavior relies on the study of periodic solution of the free system; this can be 
S^ . considered as an extension of the use of normal modes of a linear free system in order 

to study the dynamics of the forced associated linear system as explained in [GR93] (in 
English see [GR97J). Periodic solutions for non linear o.d.e. arising in structural dynamics 
have been considered in [BB051 IBGKMOH [JR091 IKPGV091 IJPS041 IABBG061 IMiklOj and 
in many others. 

In |JR10j . the Lindstedt-Poincarre method was used in order to derive approximate 
non linear normal modes (a periodic solution close to a linear normal mode) for small non 
linearity. 

At this point it is worth to recall a fundamental remark of Henri Poincarre in his report 
" Sur le probleme des trois corps et les equations de la dynamique" : 

"Ce qui rend ces solutions periodiques aussi precieuses, e'est qu'elle sont, 
pour ainsi dire, la seule breche par ou nous puissions penetrer dans une place 
jusqu'ici reputee inabordable" ( [Poi99| ). 

In English 
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"What renders these periodic solutions so precious is that they are, so to speak, 
the only breach through which we may try to penetrate a stronghold previously 
reputed to be impregnable" 

Although this remark is mainly devoted to celestial mechanic, it is also meaning-full for 
vibration of structures. 

In this report we study periodic solutions of non linear systems of ODE arising in 
free vibrations of structures subjected to unilateral springs; the ideas and the methods 
may be used for systems modeling other situations in particular when the non linearity 
is only Lipschitz. We recall a simple modeling process in section [2j Then in section [3l 
we start by presenting some ideas of quasi Newton type algorithm for solving equilibrium 
of structures in the static case with some Lipschitz non linearities; for example, these 
non linearities arise in cases where structures are subjected to unilateral springs; these 
unilateral springs are often simplified models of structures like bumpers usually made of 
viscoelastic materials like in [HFR09J; the same model will be addressed in the dynamic 
case. We derive a constructive proof of existence of periodic solutions for a one degree of 
freedom system in section 2J then the case of several degrees of freedom is addressed in 
section [5j 

Vibrations with non linearities have been considered from an experimental point of 
view in |DOL + 03[lvn)4llDET02llMCG02|IZFGC0. c )| : asymptotic expansions and numerical 
methods have been used in |RV04al lRV06l IRV04bl IRV051 IVR051 [RVM IRou08l IRoulOl 
[JR091 IBR091 IBTal URlOl IBralOj . 

The work of Lyapunov |Lya49| is often cited as a basis for the existence of periodic 
solutions which tends towards linear normal modes as amplitudes tend to zero; this result 
uses the hypothesis of analycity of the function involved in the differential system. As we 
address the case where the non linearity is only Lipschitz, this hypothesis of analyticity is 
obviously violated. 

Some general existence results, based on calculus of variations, for convex Hamiltonian 
systems are presented in [Ekc90] ; the case of analytic functions in connexion with the use 
of normal forms is considered in [JL91, TA06| IILlO| . Energy pumping has been addressed 
in several papers, for example in [SL10]. 

The case of vibration of structures with unilateral springs or more generally systems 
modeled with Lipschitz non linearities are considered in some recent studies (experiments, 
asymptotic expansions, numerical computations): [Hazl HR09a, [HR09bj [HR081 lHFR 09. 
IHazlOl [JR091 iJRltJl IJPS041 IVLP08J . A study of dynamics of elastic shocks is in |ACR02j . 
The case of distributed systems modeled by partial differential equations has been ad- 
dressed for example in [SP94J; general references for perturbation methods are among 
many others [SV85, Nay81| for differential equations and [KC81| for partial differential 



equations. 

A review of periodic solutions of non autonomous ordinary differential equations may 
be found in [Maw09] for details, see also [RM73]. 

In the static case, for example, a stability result concerning the obstacle problem for 
a plate is considered in [PL08J. 

Non smooth optimization is an active field of research; we only cite some references 
connected to this research: |DR09| . [Kun08| . [Iei04] and the references therein. 

Here we intend to use methods from non smooth optimization to prove existence of 
periodic solutions ans to derive a numerical algorithm to find them; the precise description 
and implementation of the algorithm will be addressed in a forthcoming paper. 



2 Equations from structural mechanics 

Consider a spring-mass system with some weak broken springs which are acting only in 
compression. Let u denote, the displacement of the masses, and jj the strain of the 
springs; it is related to the displacement by 7* = Uj+i — Uj or in vector form with an 
incident matrix B 

7 = Bu 

if we assume that the material is elastic linear, the stress-strain law is 

a = £7 + eE'(j + d)- 

where e is a small parameter; d denotes some backlash; in other words for some springs in 
small compression or in traction, there is no induced stress. 

In the static case, the force applied to each mass is denoted by Y and we have the 
equilibrium equation 

B T a = Y 



or more explicitely B T E~f + eB T E'(^ + d)_ = Y or 

B T EBu + eB T E'{Bu + d). 



Y 



(2.1) 



Remark 2.1. • Obviously this system is quite general; many situations of structural 
mechanics may be cast in this system by using finite elements. Moreover many other 
physical systems can be modeled with such non linearities. 

• To solve this system, the idea is to use some quasi Newton method taylored to 
systems involving Lipschitz functions. 

• The proof of convergence relies on a fixed point method. 

The principle of the method may be found in the book of Dontchev-Rockafellar 
[DR09] as a way of proving the existence of an inverse of a Lipschitz function / or 
an implicit theorem for Lipschitz functions. 

• Here the function f(X) is the left hand side and we want to solve f(X) = Y for 
some values of Y; note that here / is only a Lipschitz function; it is not continuously 
differentiable! 



Example with 5 masses and 4 springs on a straight line The incident matrix is 
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As we have enforced no boundary condition, there is a rigid body movement: Uj = c for 



j = 1, ... 5 gives the strain 7 = 0; for example, we can enforce uq 
attached) and remove it from the degrees of freedom; so we get 



(the mass is 
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Figure 1: 5 masses and 6 springs broken at supports 



In the case where the only broken spring is the second: 



B+E'iBu + d). 



-E' 2 (u 3 -u 2 + d)- 
E' 2 (v,3 -U2 + d)- 







B 



Example with weak broken springs at supports on a straight line We consider 
the case where the springs are broken at both ends; then for n masses, we have n + 1 
springs; 

71 = u\, 7j = Uj - Uj-i, 7„+i = -u n (2.2) 

For n = 5, we get (see figure [2]) 

10 
-110 
0-11 
0-11 
0-11 
0-1 

In the case of broken springs at both ends, 

a\ = E[(>yi + di)_, ctj = Ejjj, for j = 2. . . 5, a 6 = E' 6 (jq + d G 

so that, the equilibrium equation are: 

Ku + B T E'{Bu + d)_ =F 

with 

K = B EB, the rigidity matrix, and 

[B T E'(Bu + d)-]x = E[( Ul + di)_ 

[B T E'(Bu + d)-} j = 0, for j = 2, . . . 4, [B T E'(Bu + d)„} 5 = E' 6 (-u 5 + d 5 ). 



(2.3) 

(2.4) 

(2.5) 
(2.6) 

(2.7) 
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Remark 2.2. 1. We notice that £T is not invertible; see Remark 

2. Obviously, we can easily extend this example to a bi or tridimensional array of 
masses and springs, with unilateral springs or to non linear springs with lipschitz 
non linearity. 

3 From a classical Newton method to quasi Newton method 
with a Lipschitz function 

In this introductory section, we recall a constructive convergence proof of the Newton- 
Raphson algorithm based on the fixed point method; then we extend the approach to a 
quasi Newton method for a smooth function; finally, we extend it to a sum of a smooth 
and a lipschitz function. 

3.1 Newton and Quasi Newton methods in the smooth case 

3.1.1 Newton method 

When / is differentiable, in order to solve an equation f(x) = y, the idea is to approximate 
/ by its differential: we approach f(x k +i) by f(xk) + f'( x k)(x k +i — %k) so that x k+ \ is 
the solution of: 

f(xk) + f'(x k )(x k+1 -x k ) = y 

To prove the convergence, we consider the equation for the previous index 

f(x k -i) + f(xk-i)(xk - x k -i) = y 
and manipulating, we obtain 

f'(x k )(x k+ i - X k ) = -[/(Xfe) - /(iCfc-l) - f'(x k -l)(x k - Xk-l)] 

or 

(x fe+ i - Xk) = -[/' (a^r 1 [f(x k ) - /(xfc-i) - f'(x k -i)(x k - x k -i)] 

Now, it is clear that to show convergence of this sequence, we need mainly the 
Assumptions 

• ||/'(^) _1 || < c in a neighborhood of the solution; 

• and an easy consequence of differentiability of / in a neighborhood of x 
for any e > 0, there exists 5 such that for \\x k — x\\ < 5, \\xk+i — x\\ < 5 

\\f(x k ) - /(xjfe-i) - f'(x k -i)(x k - x fc _i)|| < e\\x k - x k -i\\ 



we choose e = %- and we obtain 

IKxk+i -x k )\\ < -||(x fc -x fe _i)|| 

• which proves convergence (if the initial data is close enough to the solution itself 
close enough to the given x). 

• In this method, it is well known that the main trouble is the inequality 

\\f\x)- l \\<c 

if c is large, we have to start very close to the solution in order to converge; the "tan- 
gent subspace" to the graph of / near the solution should not be "too horizontal" . 
A way to circumvent this trouble is to embed the Newton method in a continuation 
process. 

3.1.2 Quasi Newton method 

The idea of a quasi Newton method is to evaluate the derivative of / at a fixed point x. 
Still considering the case where / is differentiable, instead of I3.1.1|, we use 

f(xk) + f'(x)(x k+1 - Xk ) = y 

in which the derivative of / is evaluated at the fixed point x. To prove convergence, we 
consider the same equation for the previous index 

f(xk-i) + f(x)(x k - x k -i) = y 

manipulating, we get: 

f{x)(x k+1 - x k ) = -[f{x k ) - /Ofc-i) - f'{x)(x k - x k -x)] 

We use the same assumption on the inverse of the derivative /' and still with differentia- 
bility in the neighborhood of x, we get for any e > 

Wf(xk) ~ /(zjfe-i) - f'{x)(x k - ac fc _i)|| < e\\x k - x k -i\\ 

then, from 

(xfe+i - x k ) = -[f'(x)}^ 1 [f(x k ) - f(x k -x) ~ f'(x k -x)(xk ~ Xk-x)] 

with the same assumptions, we still get 

\\(Xk+l -Xk)\\ < -||(x fc -x fe _i)|| 

which proves convergence (if the initial data is close enough to the solution itself close 
enough to the given x). 



3.2 Quasi Newton method in a Lipschitz case 

We use some notions introduced in Dontchev-Rockafellar [DR09J. 

We no longer assume that / is differentiable but that there exists a strict estimator 
h of / at x: 

Definition 3.1. 1. We recall that the Lipschitz modulus of e at x is 

,. / -x ,. ||e(x)-e(a; / )|| 
lip(e,x) = nmsup n m 

x,x'^x,x^=x' W^ •*' II 

2. h is a strict estimator of f at x, when 

f{x) = h(x) + e(x) 
with the Lipschitz modulus of e satisfying 

lip(e,x) < fj, < +oo 
We assume that 

y = f{x) 

Moreover, we assume that h is invertible around y for x 
and 

lip{h' l ,y)<K (3.1) 

Remark 3.1. It should be pointed out that it stems directly from the previous definition 
that the assumption, the Lipschitz modulus at x of a function be finite, implies that this 
function is Lipschitz in a neighborhood of x. 

3.2.1 Examples 

• f(x) = (x)_, lip(f,x) = 1 for x < but lip{f,x) = for x > 

• e(x) = X- — ax with < a < 1 

lip(e, x) < max(a, 1 — a) 

• e = (Bu + d)- — A(Bu + d) where A is a diagonal matrix with < Xj < 1. 
set dj(u) = J2 BjkUk + dj 

as ej(u) — ej(u') = dj- — Xjdj, we get 

\ej(u) — ej(u')\ < max(Aj, 1 — Xj)\dj(u) — dj(u')\ < max(Aj, 1 — Aj)||-B||||u — u'\\ 

lip(e,x) < maxmax(Aj, 1 — Aj)||JB|| 
j 

consider now the case with 

f(u) = B T EBu + eB T E'{Bu + d)- . (3.2) 

It is the left hand side of our example from structural mechanics; as the first term is 
linear, the function e is e = e[(Bu + d)- — A(Bu + d)] and 

lip(e,x) < emaxmax(Aj,l — Aj)||-B|| 
j 



h(u) = B T EBu + eB T E'k(Bu + d) (3.3) 

is invertible when the rigidity matrix K e = B T EB + eB T E'B is invertible and its 
inverse is: 

h-\y) = K- 1 [y-B T E'Ad] 

and 

lip{h-\y)<\\K- l \\ 

So when K e is invertible, for e small enough, 

lip(e, x)lip(h~ ,y) < 1 

this property is crucial for convergence but it should be pointed out that Kq = B 1 EB is 
not required to be invertible. 

3.2.2 Algorithm 

Now the idea is to build a sequence to solve 

f(x) = V 

by using the strict estimator h: 

h(xk+i) + e{x k+ i) = y 

to compute x^+i solution of this equation is as difficult as the initial problem but it may 
be approximated by replacing e{xk+i) by e(xk) 

h(x k+ i) + e(x k ) =y 

to prove convergence, we use the formula at the previous iteration 

h(x k ) + e(xfc_i) = y (3.4) 

we get 

x k+ i = h' 1 (y - e(xfe)) x k = h' 1 (y - e(x fc _i)) 

and 

Xk+i = %k + h^(y - e(xfe)) - hT x {y - e(x k -i)) 



• with the assumption (|3.ip . we get that for e(x k ) and e(x k ^i) small enough and y 
close enough to y, 

||^ _1 b - e(x fc )] - h~ l [y - e(x fc _i)] < A|| - e(x k ) + e(x fc _i)|| 

• on the other hand, h being a strict estimator: 

\\e(x k ) - e(xfc_i)|| < u\\xk - x k -i\\ 

• moreover, we assume that [i k < 1 so that for x k and x^_i close enough to x, we 
may assume that A v < 1 



• From [3X21 we get 

||xfc+i - x k \\ < A u\\xk - xt-i\\ 

Usually, we have to show that if xq is close enough to x, x\ lies in a ball in which 
h^ 1 is Lipschitz; here it is obvious; it proves the geometric convergence when 

Az; < 1 

We have obtained 
Proposition 3.1. Consider the equation 

f{x) = y (3.5) 

Assume that 

f(x) = h(x) + e(x) 

with the Lipschitz modulus satisfying 

lip(e,x) < fi < +oo lip[K~ ,y) < k < +oo 

with fiK < 1, f(x) = y and e(x) = 
then the sequence 

h(x k+1 ) + e(x k ) =y 

with xq is given close to the solution, converges geometrically to the solution of (]3.5p close 
to Xo. 

This is the case for our example from structural mechanics when e is small enough. 

Continuation process This method may be embedded in a continuation process: 

• Set a step size 5e and an initial value of e\ and a starting point x\\ 

• For 1 from 1 to L 

• Use the non smooth quasi Newton method to find the solution x\ of the equation 
for ei with the starting point x/ ; 

• increase ei + \ = ei + 5e and set x/ + i = x*[ 

• end For 

Extensions for differential equations are considered below; control problems are possible 
extensions. 

An example Consider the iterations (|3,4h for the example of broken springs at both 
ends ([2]), we get with ([32D and (l3~3|) : 



Ku k+1 + eB T E'K{Bu k+1 + d) = -e [B T E'{Bu k + d)_ - B T E' K{Bu k + d))) + F (3.6) 

We notice that the matrix of the left hand side is the sum of K and eB T E'AB; it is usually 
invertible. 



Remark 3.2. Had we use the "natural" iterative method 

Ku k+1 = -eB T E'A(Bu k + d) + F (3.7) 

we notice that the matrix K is not invertible so that this algorithm is not working at 
all; in cases where the smallest eigenvalue of K is not zero but small as in problem arising 
in large aerospace structures, the algorithm will converge for e quite smaller than with the 
previous algorithm fj3.6f) 

It could be proved that the matrix is better conditioned with the "quasi Newton" 
algorithm than with the natural algorithm. 

4 Periodic solution of a one degree of freedom spring-mass 
system with Lipschitz non linearity 

4.1 The differential equation, periodic solution 

We set 

~ dx 

X = Tt (41) 

and we consider a one d.o.f. spring-mass system with a linear and a weak linear unilateral 

spring: 

x — \x\ 
mx + kx + ek(x)- = with x_ = (4-2) 

or more generally a linear and a weak non linear spring 

mx + kx + ekg(x) = (4.3) 

We set to 2 = — , so this equation may be written 

x + u 2 x + eu] 2 (x)- = (4.4) 

As we suspect that the frequency of a periodic solution depends on e, we perform the 
change of variable: 

6 = LO e t 

with 

1 1 - erj(e) 



UJ 2 CO 2 

Remark 4.1. This is a classical change of variable, e.g. see |Ver90| ; it is used in the 
method of strained coordinates also called Linsted-Poincarre; in particular, it has been 
used in [JR10] to derive an asymptotic expansion for a quite similar system of differential 
equations. Instead, here we use it to derive a constructive proof of existence of a periodic 
solution; we are only looking for periodic solutions whereas in [JR10] , we can obtain a 
quasi periodic expansion. 

and we set 

If we set x(9) = x(t), equation (|4.4p may be written: 

x" + (1 — er])x + e(l — er/)x_ = and in the general case (4-6) 

x" + (1 - erj)x + e(l - erf)g{x) = (4.7) 

10 



or 

x" +x + ef(x,i],e) = (4.8) 

with 

f{x, 77, e) = — rjx + (1 — €7])x- and in the general case (4-9) 

f(x,Tj,e) = -r)x+(l-er))g(x) (4.10) 

As the solution of (14.81) is: 



x = acos(9) + bsin{6) — e / sin(9 — s)f(x(s),r],e)ds (4-11) 

or 

x = acos{6) + bsin(9) - er(6, a, e, 77, b) ( 4 -12) 

we obtain the following lemma which may be considered as an extension to the Lipschitz 
case of a classical result which, for example, may be found in [Ver90j. 



Lemma 4.1. Consider the solution x of the Cauchy problem of equation (|4.2|) (resp. 
) and the solution x of the associated equation (|4.8j) . (|4.9p (resp. (|4.10p ) after change of 
variable. For e close to zero, x is a solution of (|4.2|) (resp. (|4.3p ) of period 2tt/uj € if and 
only if x is a solution of ()4.8|) . (|4.9j) (resp. (|4.10p ) and 

r(27r,a,e,r],b) = and ' — - — =0 or equivalently F(p,y) = (4.13) 

with 

p= (a,e), y = (rj,b) (4.14) 

r-2ir 

Fi(p,y)= sin(s)f(x,r],e)ds (4.15) 

Jo 

i-2-it 

F2{p,y) = I cos(s)f(x,rj,e)ds (4.16) 

Jo 

4.2 Computation of 77(0) 

Case x- We have 

/(x, 77, 0) = — r/x + X- 

and for 6 = 0, 

f(x,rj,0) = —nacos(6) + (acos(#))_ 

so 

f-2n 

F 1 = / sm(s)(acos(s)))_ds (4.17) 

F2 = —rjacos(s) 2 + cos(s)(acos(s)))-ds (4-18) 

Jo 

we note that F\ is identically zero; we remark that 

if a > 0, (acos(s))_ = a(cos(s))_ and if a < 0, (acos(s))_ = a(cos(s)) + (4.19) 
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in both cases, we get from F 2 

r?(0) = ~ 
by using for example 

2 4 + °° (— T) fe 

I cos ( s )l = z ~ z E i^3i cos ( 2 ^)' ( 4 - 2 °) 

fc=i 

and by anticipating that 77 is a Lipschitz function of e (see proposition 14.21 ) 

J_ _ 1 - e/2 + o(e) 

w 2 w 2 

and 

x = 0005(6!) + bsin{9) - er(9, a, 0, ry(0), 0) + o(e) 

This result might be obtained by direct inspection of the level curves of the energy as- 
sociated to the differential equation, see [JRlOj; but this approach may be extended to 
systems of differential equations. 

General case For simplicity, we consider b = 0, the equation F\{p,y) = is identically 
satisfied. The second equation yields: 

1 f 27V 
77(0) = — / cos(s)g(acos(s))ds (4.21) 

4.3 Perturbation of the solution x with respect to y — (77, b) 

Lemma 4.2. Assume that x is solution of (|4.8p (I4.10p with f and g Lipschitz with respect 
to all variables; 

\g{x 2 ) - g{xi)\ < k\x 2 - xi\ (4.22) 

1/(^2, 772,62) ~ /(^l,^l>ei)l < k{\x 2 -x\\ + I772-771I + |e 2 -ei|) (4.23) 

assume that the initial data are: 

x(0) = a a , x'(0) = eb a , (resp. x'(Q) = b a ) with a = 1,2 (4.24) 

then, x and F are Lipschitz with respect to y = (77, b) and p with a modulus of magnitude 
ek{resp. k) (where k < +00 ): 

V0€[O,2tt], \\x 2 (6) - xi(6)\\ < ek(\b 2 - h\ + \7]2 - m\) and (4.25) 

resp. when the initial velocity is not of order e (4.26) 

V#€ [0,2vr], \\x 2 {6) - x x {9)\\ < k(\b 2 -h\ + \r]2-m\)], and (4.27) 

V^€[0,2vr], \g(x 2 (6)) - g( Xl {e))\ < ek(\b 2 - h\ + \m ~ m\) (4-28) 

but we have only (4.29) 

\f(x(e,y 2 ), m ,e) - f(x(6, yi ), m ,e)\ <k{e\b 2 - &i| + \n 2 - h Y \) (4.30) 

\\F(p,y 2 ) ~ F{p, Vl )\\ < k{\b 2 - h\ + |t/ 2 - 771I) (4.31) 
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Proof. The proof relies on the formula (|4.1ip . we get: 



x 2 (0) - X!(0) = e(6 2 - bjsinid) - e / sin{6 - s) [f(x 2 , r? 2 , e) - f(x 1: m , e)] ds (4.32) 

JO 

|x 2 (0) - xi (0)1 < e [|6 2 - fei[ + 2vrA:|r/2 - ??i|] + 2vreA;|a;2 - x t \ (4.33) 

from which we get equation of (|4.25|1 and we deduce from (|4.23p 

\g(x(9,y 2 )-g(e,yi)\ < ek(\b 2 - h\ + \m~m\) (4-34) 

which is (|4.28p ;but we get only 

\f(x(e,y 2 ), m ,e) - f(x(e,yi), m ,e)\ <k(e\b 2 - h\ + \rj 2 - h x \) (4.35) 

we deduce equations (|4.3ip . D 

4.4 A fixed point approximation method 

Consider now the solution of (|4.3p , with / given by (|4.10p ; when g(acos(s)) is an even 
function, a 2n periodic solution with 6 = satisfies trivially F\ = 0; with lemma |4T2| we 
obtain: 

Proposition 4.1. For a a 2n periodic solution with 6 = satisfies 

r\ = (1 — er])F(a, e,r?) (4.36) 

with 

TV \ fo n cos(s)g(x(s))ds 

J cos(s)x(s)as 
and for e small enough, this function satisfies 

\F(a,e,7] 2 ) -F(a,e,rn) < ek\r] 2 - rji\ (4.38) 

so that the following sequence 

v k+1 = {l-er) k )F(a,e,r) k ) (4.39) 

converges and it proves the existence of a 2ir periodic solution of (|4.8p with (|4.10p . hence 
the existence of a periodic solution of (|4.3p with angular frequency (|4.4p . 

4.5 Estimator problem 

We consider another approach which should be better conditioned (see Remark 13.21 ) and 
enable to consider arbitrary initial velocity whereas in the previous paragraph, the initial 
velocity is zero. 

Equations (|4.13p may be considered as defining an implicit function 

p t— > y (4.40) 

As in previous subsection, this function F is not smooth and we are going to show that 
the implicit function (|4.40p may be defined with a fixed point of a contraction maping 
following general lines of [DR09]. So we introduce a differential equation and a function 
H that will be proved to be an estimator of F. 
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Definition 4.1. Consider the differential equation 

e , +Z + eh(Z, V ,e) = (4.41) 

with 

h(x,r],e) = -r)£+ (1 - eifja^ with < a < 1 (4.42) 

and we introduce the following function: 

H 1 (p,y)= sin(s)h(£,ri,e)ds (4.43) 

Jo 

[■2lT 

H 2 (p,y)= cos(s)h(£,r],e)ds (4.44) 

Jo 

note that for e = 0, h = f so that H = F. 

Method: to solve 

F(p,y)=0 

following the method introduced in paragraph 13. 2\ we use the sequence y k , defined by 






H(p,y k+L ) = -E(p,y k ), where E(p,y) = F(p,y) - H(p,y) (4.45) 

As H is a smooth map, we prove that H is invertible for e small enough by using the 
classical implicit function theorem (see e.g. . 

Lemma 4.3. 1. The Jacobi matrix of 

y — > H(p,y) 
with p = (a, e) and y = (r],b) at e = and b = is 

(-77(0) + a)7r 

-an 

2. so with a 7^ 77(0) = | and /or e small enough, the equation 

H(p,y)+e = (4.46) 

/ias a solution y(p, e) which is differentiate and so locally Lipschitz with respect to 
its variables p,e for e small and p close to (a, 0). 

Proof. For e = and 6 = 0, we notice that 

^ = m = sin(e) 

d V ' db [ h 

so 

9h ,„, dh 

— = -acos(V), — = {-n + a)sin(9), 

from which we deduce the first part of the lemma. The second part is deduced from the 
classical implicit function theorem. □ 
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Remark 4.2. With examples 13.2.11 when / is defined by (|4.9p and h by (|4.42p then f — h 
satisfies the general inequality (|4.4T[) . 

Lemma 4.4. Assume that x is solution of f|4.8j) (|4.10jl and £ solution of (|4.4ip . (J4.42J) 
with the same initial conditions and with f and g Lipschitz with respect to all variables 
and moreover that h is a strict estimator of f , i.e.: 

For small £, |/(£,T/,e) - h^,r,,e)\ < /z|£| (4.47) 

i/ien 



|x(0)-£(0)|<e/xc Sup (|£|) (4.48) 

0<s<2tt 



/or € [0, 2tt] 



Proof. Indeed, from (|4.11j) with x and a similar formula for £, both having the same initial 
conditions, we get 



x (0) - £(0) = e / sin(0 - s) [/(a?, 77, e) - /i(f , r/, e)]ds (4.49) 

Jo 

and so (4.50) 

x{6) - m = e [ sin(6 - s) [f(x, r), e)_/(£, 7?, e) + /(£, r,, e) - /»(£, r), e)] (4.51) 

JO 

with (IQTjl and / Lipschitz, |ar(0) - f (0)| < efc / (|s - CI + ju|f |)d« (4.52) 

JO 

with Gronwall lemma (jl2), we get 

\x(9) - £(0)| < enSup(\Z\)e 27rek (4.53) 

for e [0, 2vr] D 

Lemma 4.5. Assume (|4.47j) . i/ien i/ie function E(p,y)=F(p,y)-H(p,y), satisfies 

lip(E, y) < e u < +oo (4.54) 

wi/j po = ( a i0) T wii/i a arbitrary, in other words, H is a strict estimator of F, uniformly 
in p for e < e (following "DR09^) . 

Proof. To prove it, we use the short hand / = f(x(9,y),fj,e), h = /i(£(0,y),7], e) where 
y = (fj,b) is close to y = (77, b). We have to estimate 

E(p, y) - E(p, y) = F(p, y) - F(p, y) - (H(p, y) - H(p, y)) (4.55) 

where F is defined in (|4.15p . (|4.16p and H in ()4.43j) . (j4.44|) . so the basic point is to consider: 

/ - / - (h - h) = - r)x + (1 - eff)g(x) 

- (-r)x + (1 -677)51(0;)) 

(4.56) 



- {-VC + (1 - efj)g(0 
+ (-r£ + (l-eri)g(Z)) 
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Then we split the right hand side into 3 expressions that we manipulate separately 

\v(x - - fj(x - |)| =\ti(x -£-£ + £) + (t}- fj)(x - £)\ 
<rj[\x - x\ + |f — £| + \x - £|| 7/ — fj\] 
with lemma I4T21 and |4~31 

<ec(l + Sup\£\)\y - y\ 



Similarly 



and 



\g(x) - g{x) - g(£) + g(0 \<k\\x-x\ + \£-£\ 

< ek\y-y\ 



(4.57) 



(4.58) 



-fjg(x) + rjg(x) + f]g(£) - rjg(£,) 

iv - fj)g(%) + v(ai x ) - 9{x)) - (v - v)g(Q - fj(g(Q - g(€) 



< € 



{Sup\g(x)\ + Swp\g{£)\)\q -f}\ + Sup\fj\(\g(x) - g{x)\ + \g{£) - g{£)\) 



< ec\y-y\ 

The last inequality is a consequence of lemma 14.21 and | 
So we have obtained 

\f-f-(h-h)\<en\y-y\ 



(4.59) 



(4.60) 



uniformly with respect to p = (a, e) from which the lemma is obtained using the definitions 
of E(p, y) and E(p, y) where F is defined in (j4~T5"j) . (liTTHjl and H in P~l3]) . p~4"4]l . D 



Proposition 4.2. Equation (J4.8J) . (I4.9P (resp. (I4.10p J /ms a 27r periodic solution for e 
small enough and it may be computed with the iterative process (J4.45J) ; moreover, rj,b are 
Lipschitz function of a, e and so 

\x - acos(6) - bsin(9) + er(9,a,0,n(0),0)\ < ke 2 (4.61) 

Proof. The proof is now simple: write the iterative method (J4.45P for k and k — 1 



H( P ,y k+1 ) 
H( P ,y k ) = 



---E(p,y h ) 
-E(p,y k - r ) 



and by subtraction 

H(p, y k+1 ) - H(p, y k ) = E(p, y^ 1 ) - E(p, y k ) 
with previous lemma I4T31 there exists A 

ll/^-Zll^AH^p,/- 1 )-^,/)!! 

and as E is a strict estimator (lemma I4.5P 

||/+i_/||< eA /,||/-/- 1 || 



(4.62) 
(4.63) 



(4.64) 



(4.65) 



(4.66) 
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As for the static case we should prove that for p close to po and y° close to the solution 
obtained for e = 0, y 1 lies in a ball where y(p,e) defined in lemma I4T31 is Lipschitz. 

The first part of the proposition is then proved by the fixed point theorem. 

The fact that y is a Lipschitz function of p is derived with a somewhat classical argu- 
ment. With lemma H31 we get from 

H{p 2 ,y 2 ) + E(p 2 ,y 2 ) = 0, H( Pl , yi ) + E( Pl , yi ) = (4.67) 

\y2-yi\<\\E(p 2 ,y2)-E(p 1 ,y 1 )\ (4.68) 



and with lemma [4"31 

\V2 -yi\ < K k \P2 ~Pi\ + e/%2 -yi\) (4.69) 

from which 

\y2-yi\<- (fc|p2-pi|) (4.70) 

1 — efj, 

Then, as r is Lipschitz with respect to p and y, we obtain the inequality (|4.6ip D 

5 Vibrations with several degrees of freedom 

5.1 Differential equation model of a spring-mass system 

In the dynamic case, equation (|2.ip becomes: 

Mu + B T EBu + eB T E'(Bu + d)-=0 (5.1) 

with 



(5.2) 



du 

u = — - 
dt 

Remark 5.1. Many other mechanical systems with unilateral properties can be cast into 
this frame in particular when using finite elements see e.g. |JPS04] ; the case of beams with 
support on a unilateral spring is considered in H Hazim |HazlO| . As for one degree of 
freedom, we can consider a more general case: 

Mil + B T EBu + eQ(u) = (5.3) 

with Q(u) a Lipschitz function; as an introduction we shall also consider the case of G(u) 
a smooth function; it can model a non linear stress-strain law. 

5.2 Eigenvectors, estimator system 

In order to derive approximate periodic solutions, it is convenient to write the differential 
system in an eigenvector basis; we denote K = B EB the rigidity matrix and introduce 
the matrix of generalized eigenvectors $ = [$i, . . . , $j, . . . , <& n ] where <!>/% are generalized 
eigenvectors associated to generalized eigenvalues io\ satisfying 

- wgM$ fc + K$ k = (5.4) 

or in matrix form: 

K§ = n 2 M$ 
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with Q, 2 the diagonal matrix of eigenvalues ; we assume 

$ T M$ = I 

consider the change of function 

u = <3?x 

we obtain: 

x + tt 2 x + e$ T B T E'{B<&x + d)_ = (5.5) 

and in the general case 

x + n 2 x + e$ T £(<M) = (5.6) 

For unilateral spring, as 7 1 — > f('j) = 7- is a concave, Lipschitz function, we can 
use as strict estimator (in the sense of [DR09J) 7 1 — > h{^) = A7 with A = diag(Xj), 
< Aj < 1; if we set e = / — h, we have: 

lip(e,0) = maxmax(Aj,Aj — 1) 
3 

so a reasonable estimator differential system is 

C + n 2 £ + e$> T B T E'k{B^ + d) = (5.7) 

5.3 Periodic solutions 

In the linear case when e = 0, equation (|5.ip has periodic solutions of the form: 

«fc = (« fc e* Jk *+15fce- <w **)* fc (5.8) 

where 3>& are generalized eigenvectors defined in (15. 4p . 

5.3.1 Orientation 

We intend to prove that for e small enough, there exists periodic solutions to equation 
(|5.ip . or in the general case (|5.3p : the method is constructive and will enable to derive a 
numerical scheme to approximate it; we use the equivalent system (|5.5p or in the general 
case (|5.6|) and we introduce the following notations: 

p = (ai,e) T y = (r],b 1 ,a 2 ,b 2 ,...,b n ) T (5.9) 

Remark 5.2. Notice that a\ plays a particular role: we are going to consider a periodic 
solution close to the harmonic solution of period -^- of (|5.5p with e = 0; obviously, any 
other harmonic solution could be chosen by using a permutation of indexes. 

5.3.2 For unilateral springs 

Following the lines of the 1 d.o.f. case, we introduce a change of variable suited to look 
for periodic solutions close to the eigenmode of index 1. We set 9 = uj € t with 



1 1 — eri(e) , dx 

-^ = 5 — and set x = — 

uit 0Ji A dv 



x(6) = x(t), — = ]±-L and set x' = — (5.10) 
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Equation (|5.5p becomes 

x"i +xi- er]xi + e ~ ^ ^T B T E' (B^x + d)_ = (5.11) 

xy" + ^Xj - flj^ay + e ~2 V <i>JB T E'(B$x + d)_ = 0, j = 2, . . . , n (5.12) 

U^ Wj (J-^ 

for the estimator differential system, 

£i" + £i - er& + e S}-ip±<S>T B T E§1 B T A{B<S>£ + d) = (5.13) 



a; 



,2 , ,2 



- , ^e ^3c , e ( 1 " er ?)^T D T 



6" + 4£i " «/4^i + V 2 " $j £ J £A(.B$£ + d) = 0, j = 2, . . . , n (5.14) 



WJ wj w 1 



we set: 



fj x , n, e) = -nx! + ^-^<5>?B T E'(B<S>x + d)_ (5.15) 

fdx,r],e) = -r!%x j + ^^-^B T E'{B^x + d)^ = 0, j = 2,...,n (5.16) 

so that the system may be written: 

x"i+xi + e/i(a;,?7,e)=0 (5.17) 

x "i + ^2 x j + ^fj{x,il,e) = 0, j = 2, ...,n (5.18) 



similarly, we set: 



hi(S,ri,t) - -^i + 1—P-$Tb t E'A(B$Z + d) (5.19) 



o,'T 



,2 



/i i (e,r ? ,e) = -r ? ^i< j + ^^^jB T ^ / A(S^ + d) = 0, j = 2,...,n (5.20) 

and the "estimator" system is now: 

6"+6 + eM^,e) = (5.21) 

G" + 4&+ e M£,^)=0, J=2,...,n (5.22) 

5.3.3 General case 

In the general case of equation (|5.3[) , in the eigenvector basis (15 .6p , we still obtain equation 
(l5T?l) . ([57L8t but with 

h(x, rj, e) = -rjx! + ^—P-®1 G($x) (5.23) 

CO 2 1 — 

/ i (x,r ? ,e) = -r ? 4a; J + £LQJg($ x ) = 0, j = 2,...,n (5.24) 

in the non smooth case, we consider a strict estimator /ij(x) of QjQ(Qx) and the estimator 
system is written as previously in (I5.2ip . (l5.22p 
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5.3.4 Condition of periodicity 

Lemma 5.1. Consider the solution x of the Cauchy problem of system (|5.5[) or (|5.6|) and 
the solution x of the associated system (|5.17p . (J5.18P after change of variable and with 
initial conditions: 

x 1 (0)-a 1 = 0, xi(0)-6i = (5.25) 

: : (5.26) 

x n (0) - a n = 0, x n (0) -b n = (5.27) 

(5.28) 

For e close to zero, x is a solution of (|5.5p of period — with —% = — ^^ if and only if 

F(p,y) = (5.29) 

where (p, y) is defined in (|5.9p and the function 

F:I 2 xR 2 "4 R 2n (5.30) 

is defined by: 

f-2n 

Fi(p,y)= sin(s)fi(x(s),rj,e)ds (5.31) 

Jo 

/>2tt 

F2(p,y)= cos(s)fx(x(s),rj,e)ds (5.32) 

Jo 

F 3 (p,y)=x 2 (2ir)-a 2 (5.33) 

! (5.34) 

•^2*n-i(p, 2/) =Xn(27r) - a n (5.35) 

F2*n(p,y) =x n (2-K) - b n (5.36) 

(5.37) 



Proof. The proof is similar to the proof of lemma 14.11 



□ 



In the non smooth case, it will be proved that the following function H is a strict 
estimator of F. 



Definition 5.1. Consider the solution £ of the Cauchy problem of system (|5.7p in the 
form (I5.2ip . (|5.22p with initial conditions: 

£i(0)-ai = 0, 

£n(0) - On = 0, 



£i(0)-6i = 


(5.38) 




(5.39) 


4(0) " ^n = 


(5.40) 




(5.41) 
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the function H is defined by: 

(■2-K 

Hi(p,y)= sin(s)hi(£(s),r],e)ds (5.42) 

Jo 

(■2-K 

H 2 (j),y) = cos(s)hi(£(s),rj,e)ds (5.43) 

Jo 

H 3 {p,y) =U2n) - a 2 (5.44) 

: (5.45) 

H 2 * n -i(p, y) =^„(27r) - a n (5.46) 

H 2 *n(p,y)=L(^)-b n (5.47) 

5.4 The smooth case 

As indicated in the introduction, the analytic case is well known; here we just assume 
some differentiability; this type of result may be derived by other methods;we emphasize 
that this approach is constructive in particular if we approximate x with a Fourier series 
(harmonic balance or spectral method). 

Lemma 5.2. In the case where Q is smooth, the Jacobi matrix of 

y^F(p,y) (5.48) 

is invertible at po = (a, 0) with b\ = 0; where F is defined by (|5.3ip . (|5.37p 

Proof. We get for e = 

xi = aicos(9) + bisin(9) (5.49) 

Xi = a jC os(^-6) + bi—sini^-e) (5.50) 

_ = , s =»(.| (5.51) 

^-0. |l = C STO (5.52) 



so that for b\ = 0, 



-A = - ai cos(9) -£ = -r]sin(9) + 4r*?'G / ($aO$i«*n(0) (5.53) 

or] obi u( 



so that 



dF i n dF i l xxm 

or] obi ^>{ 



-i = -rjTT + -j$f G'($x)$i 7T (5.54) 



-5— = -ai?r — - = (5.55) 

on obi 

ia ^w§*° (5 - 56) 
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For other indexes 



^l = cos(^e), ^ = — sm(^0) (5.57) 

Odj 0J\ OOj LOj U)\ 

d ll = -^esin^e), d -^ = cos(^e) (5.58) 

Odj (jj\ uj\ obj uj\ 



9F ^- dx \2,)-l, °*t± = %L(2K) (5.59) 



dcij dcij dbj dbj 

dcij dcij ' dbj dbj 



(5.60) 



dF 2j-l to U i\ 1 dF 2j-l Wi . . Uj 

— ^ — = cos(27r— ) — 1 — -7 — = — sinizn—) (5.61) 

dcij uj\ dbj ujj oji 

^k = _% m(27r ^) ^i = co S (2n^) - 1 (5.62) 

dcij loi loi dbj (jJ\ 

When -^- is not an integer for j > 1, 

the determinant 

d d( F 2*j-i,F 2 ,j) 

d6t d(a'j + \bf^) + ( } 

ans so the Jacobi of F is invertible. □ 

Proposition 5.1. Equation (I5.3ip to (I5.37P /ias a solution for e small enough and this 
solution may be computed with an iterative method. 

Proof. We conclude the first part from the classical implicit function theorem; the solution 
may then be approximated by a Newton method; it may be embedded in a continuation 
process. □ 

5.5 Strict estimator of F 

Lemma 5.3. The function function E(p,y)=F(p,y)-H(p,y), satisfies 

lip(E, y) <e/j. < +00 (5.64) 

at po = (ai,0) T with a\ arbitrary, (following the general lines of WR 091) in other words, 
H is a strict estimator of F uniformly in p close to po . 

5.5.1 Proof 

The solution of (|5.5[) after change of variable 8 = u € t, is solution of (|5.1ip and (I5.12P or 
(J5TTT1) . (J57L8]) ; it may be written 

f e 
xi =a\cos{9) + b\sin(8) — e / sin((9 — s))fi(x(s),r],e) ds (5.65) 

Jo 

[•9 

Xj =ajcos{—0) + bj — sin(—0) — e — / sin{—(6 — s))fj(x{s),n,e) ds (5.66) 

UJ\ UJj UJ\ IjJj Jq uj\ 
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and the solution of (|5.7p after change of variable is solution of (|5.2ip and (|5.22j) with 
h defined in (f5TT9l) . ([5T20]) 



£l = CLjCOS 



f 9 

(6) + bjsin(9)-e sin(9 - s)hi(£(s),r),e) ds (5.67) 

Jo 



& = Qjcosi^-O) + bj—sin^O) -e— sin(^-{9 - s))hd£{s),ri,e) ds (5.68) 
and so 

f-2n 

Ei(p,y)= sin(s)[fi(x(s),r],e)-hi(£(s),ri,e)] ds (5.69) 

Jo 

/>2vr 

E 2 (p,y)= cos(s)[f 1 (x(s),7 ] ,e)-h 1 (C(s),v,e)} ds (5.70) 

Jo 

r-T 

Ejfay) =Xj(T) - Cj(T) = -c^i / M(^(T- S ))[/ i (x( S ),» ) ,e)-/ tJ K( S ),»),e)] ds 

(5.71) 

As the solution of a system of differential equations is Lipschitz with respect to the 
initial conditions and with respect to time, we obtain the lemma from slight manipulations 
as in the proof of lemma B3| absolute values are replaced by norms and we use the following 
lemmas. 

5.5.2 Some lemmas 

Lemma 5.4. Assume that x is solution of (J5.17P (J5.18P with f and g Lipschitz with respect 
to all variables; 

\\g(x) — g(x)\\ < k\\x — x\\ (5-72) 

\\f(x,fj,e) - f{x,r],e)\\ < k(\\x-x\\ + \fj-T]\ + \e - e|) (5.73) 

then, x and F are Lipschitz with respect to y = (rj, b) and p; more precisely assume 
that the initial data are: 



then: 



x(0) = a, x'(0) = b, (resp. x(0) = a, x'(0) = b) (5.74) 

V#€[0,2tt], \\x(9) - x(9)\\ < k(\\b - b\\ + \fj - r]\) and (5.75) 

V0€[O,27r], \\g(x(8)) - g(x(6))\\ < k(\\b - b\\ + \fi - V \) (5.76) 

and we have (5-77) 

\\f(x(9, y),rj, e) - f(x(9, y), v , e)|| < k (\\b - b\\ + \rj-r,\) (5.78) 

\\F(p, y) - F(p, y)\\ < k(\\b - b\\ + \fj - r,\) (5.79) 



Proof. The proof relies on the formula (15.65P , we get: 

Xl {9) - Xl {9) = (h - b l )sin{9) - e i sin{8 - s) [h(x, fj, e) - h(x, V , e)] ds (5.80) 



Jo 
\xi(6) -xi (0)|| < [\\bi - 6i|| +27reA;|^- 7/|] + 2vre/c Sup \\x(s) - x(s)\\ (5.81) 

0<s<2tt 
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and from (|5.66p . we get 
Xj(9) — Xj(8) = — (bj — bj)sin(-O) — e / sin(— (9 — s)) [f(x, 77, e) — f(x, 77, e)] ds 

LOj L • LOi Jq LJi 



|*#)-*#)||<- 



\bj — bj\\ + 2irej\f} — r]\ + 2nej Sup \\x(s) — x(s) 

0<s<2tt 



(5.82) 
(5.83) 



from which we get equation of (|5.75p and we deduce from (|5.T3|) 

\\g(x(9,y) - g(9,y)\\ < k(\\b - b\\ + \fj- r/|) (5.84) 

which is (|5.76|) : we get 

\\f(x(6,y),fj,e) - f(x(9,y), V ,e)\ < k (\\b - b\\ + \fj - h\) (5.85) 

we deduce equations (|5.79p . D 

Lemma 5.5. Assume that x is solution of (J5.17P (J5.18P and £ solution of (|5.2ip . ()5.22|) . 
with the same initial conditions and with f and g Lipschitz with respect to all variables 
and moreover that h is a strict estimator of f , i.e.: 

Forsmallt, ||/(£,T/,e) - h(Z, V ,e)\\ < ^||£|| (5.86) 

then 

\\x(9)-m\\<WC Sup (U\\) (5.87) 

0<s<2tt 

for 9 g [0, 2vr] 

Proof. The proof is quite similar to the 1 d.o.f. case of lemma 1*01 □ 

5.5.3 Implicit equation with the strict estimator H 

In order to compute and to prove the existence of periodic solutions, we consider a con- 
structive approach as suggested in the general lines of |DR09] . 

5.5.4 Iterative method 

Solve 

H(p, y k+1 ) = -E(p, y k ) with p = [a x , e] T , (5.88) 

y k ^ = [r, k +\b^\a k 2 + \...,b k n +1 f (5-89) 

In other words, set x k (resp. £ k ) the solution of (I5T7J) ([BTTgll (resp. (|S3T]) . ([5321 ) for 
the value y = y k ; similarly set f k = f(x k ,i] k ,e), h k = h(£ k ,ri k ,e); we have to solve (15. 7|) 
for e small enough and with initial conditions (|5.38p where the initial value a\ is prescribed 
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and others are to be found as well as the frequency parameter r] such that: 

(5.90) 

(5.91) 

(5.92) 

(5.93) 
(5.94) 
(5.95) 

As H is a strict estimator of F, we may assume that for e small enough, the right hand 
side is as small as needed; we denote it a, (3: 

a\= l sin{s){f\-h\)ds, $\ = f cos(s)(tf - h\)ds, (5.96) 

Jo Jo 



i-2-k 

/ sin(s)h k+1 = 
Jo 


i-2-k 

- / sin(s)(f k - h k )ds 
Jo 


j-2k 

/ cos{s)h k+1 = 
Jo 


- ^ cos(s)(f k - h\)ds 
Jo 


e 2 +1 (27r)-a k 2 +1 = 


-[x k 2 (2ir)-$(2n)] 


e n +1 (2ir)-a k n +1 = 


-[z*(27r)-£(27T)] 


i k n + \2K)-b k n +l = 


-[%(2k) - &(2k)] 



o 



j 



x k (2n)-$(2n), /3| = x|(27r)-e|(27r) (5.97) 



Equations (|5.9Up to (|5.95p with small right hand side denoted a, j3 are 



2tt j-2-K 



sinfsjftf 1 = -ai, / cos(s)^ +1 = -/3 2 , (5.98) 

o Jo 

^ +1 (2vr) - a k+1 = -a k , ^ +1 (2vr) - b k+1 = -b) (5.99) 

Lemma 5.6. When the angular frequencies satisfy the property 

is not an integer for j > 1, 
the jacobian matrix of 

y ^H(p,y) (5.100) 

defined in (|5.42p to (J5.47P is invertible at po = (a, 0) with b\ = 
Proof. We get for e = 

£i = aicos(6>) + b lS in{6) (5.101) 

(5.102) 



(5.103) 

(5.104) 
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UJj UJl 




drj 


— = „„(*) 




drj 


— = co S (0) 



so that 



-^i = -a lC os(0) - 7 ^- = -r]sin(9) + ^^B T E'AB^[sin(e),0,...,0] T (5.105) 
a?7 ob\ 



UJ 



25 



so that 



orj ob\ uf 



— - = -r]TT + -^<&l B T E'AB$[1, 0...0] T vr (5.106) 



*<M + „ (5 ,„ 8) 

because for a suitable choice of A, we have Sp* ^ ; the other derivatives are zero. For 
other indexes 

^L = cos (£l e)i d A = ^l sin (^e) (5.109) 

Odj UJ\ OOj UJj Ul\ 

?k = -^L esin{ ^Le), d J± = cos(^6) (5.110) 

ocij uj\ (jj\ obj u)\ 

9aj <9oj ' 96j dbj 

dH 2 j-i , ujj dH 2 j-i ui\ . . loj ,,„„, 

— — - — = cos(2tt—) — 1 — — - — = — sin(2TT—) (5.113) 

ocij uii obj ujj uj\ 

dHoi uj,\ , ujj . dHo4 , ujj . 

-^-sin(2n- L ) --^- = cos(2ir^-) - 1 (5.114) 



dcij u)\ u)\ dbj u)\ ' 

the determinant 



When -2- is not an integer for j > 1 



det K ZV'::?' + (5.115) 



L d(H2*j-i,H 2*j ) 
ans so the Jacobian of H is invertible. □ 



By using the classical implicit theorem, we get 
Lemma 5.7. equation (|5.88p or (|5.90p to (|5.95p define a lipschitzian implicit function: 

p = [oi, 6] T , [a, /3] r -► y fc+1 = [T fe +\ 6*+\ a* +1 , . . . , &£ +1 ] T (5.116) 

5.6 Computation and existence of periodic solutions 

Proposition 5.2. Equation (|5,31|) to (|5.37|) has a solution for e small enough and this 
solution may be computed with the iterative method (|5.88p . 

Remark 5.3. The iterative method may be embeded in a continuation process by increasing 
the value of e 
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Proof. The proof is now simple write the iterative method f|5.88j) for k and k — 1 

H(p,y k+1 ) = -E(p,y k ) (5.117) 

H(p,y k ) = -E(p,y k - 1 ) (5.118) 

and by substraction 

H(p, y k+1 ) - H(p, y k ) = E(p, j,*" 1 ) - E(p, y k ) (5.119) 

with previous lemma, there exists A 

\\y k+1 -y k \\ < A\\E(p,y k - x ) -E(p,y k )\\ (5.120) 

and as E is a strict estimator (lemma I5.64p 

\\y k+1 -y k \\<e^\\y k -y k - l \\ (5.121) 

the proposition is proved □ 

Proposition 5.3. Equation (pU?) . (foTTB]) with fi,...f n given by ([533]) . (f5T24")) has a 2vr 
periodic solution for e small enough and it may be computed with the iterative process 
()4.45p : moreover, r\,b are lipschitzian functions of p = (oi,e) 

Proof. The proof is similar to the one of Proposition 14.21 □ 

Remark 5.4. This result should be compared with Proposition 2.1 of |JR10j : in this paper 
we have an explcit approximate value of the frequency which gives an approximation of 
the solution in an interval [0,7e -1 ]; here we get the existence of a period which gives an 
approximation of the solution on [0, oo]; but if this period is computed numerically, the 
approximation of the solution will only remain valid on a finite interval to be determined. 
The result of |JR10j may be viewed as a stability result: if the frequency is computed up 
to the order e 2 , the approximation of the solution remains valid in an interval [0, 7e -1 ]. 

6 Conclusion 

This approach seems to be the first one to provide a rigorous and constructive proof of 
existence of periodic solutions of non smooth systems of differential equations arising in 
structural mechanics. 

Implementation of the algorithm with the harmonic balance principel is in progress. 

It paves the way to reduced order modelling (see e.g. [TA06J ) for vibrating structures 
with non smooth non linearities such as contact with unilateral springs modeling bumpers. 
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